pure subroutine splder(t,n,c,k,nu,x,y,m,e,wrk,ier)
!
! calling sequence:
! call splder(t,n,c,k,nu,x,y,m,e,wrk,ier)
!
! input parameters:
! t : array,length n, which contains the position of the knots.
! n : integer, giving the total number of knots of s(x).
! c : array,length n, which contains the b-spline coefficients.
! k : integer, giving the degree of s(x).
! nu : integer, specifying the order of the derivative. 0<=nu<=k
! x : array,length m, which contains the points where the derivative of s(x) must be evaluated.
! m : integer, giving the number of points where the derivative of s(x) must be evaluated
! e : integer, if 0 the spline is extrapolated from the end spans for points not in the
! support, if 1 the spline evaluates to zero for those points, and if 2 ier is set to
! 1 and the subroutine returns.
! wrk : real array of dimension n. used as working space.
!
! output parameters:
! y : array,length m, giving the value of the derivative of s(x) at the different points.
! ier : error flag
! ier = 0 : normal return
! ier = 1 : argument out of bounds and e == 2
! ier =10 : invalid input data (see restrictions)
!
! restrictions:
! 0 <= nu <= k
! m >= 1
! t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
!
! other subroutines required: fpbspl
!
! references :
! de boor c : on calculating with b-splines, j. approximation theory 6 (1972) 50-62.
! cox m.g. : the numerical evaluation of b-splines, j. inst. maths applics 10 (1972) 134-149.
! dierckx p. : curve and surface fitting with splines, monographs on
! numerical analysis, oxford university press, 1993.
!
! author :
! p.dierckx
! dept. computer science, k.u.leuven
! celestijnenlaan 200a, b-3001 heverlee, belgium.
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
! latest update : march 1987
!
!++ pearu: 13 aug 2003
!++ - disabled cliping x values to interval [min(t),max(t)]
!++ - removed the restriction of the orderness of x values
!++ - fixed initialization of sp to real(RKIND) value
!
! ..scalar arguments..
integer, intent(in) :: n,k,nu,m,e
integer, intent(out) :: ier
! ..array arguments..
real(RKIND), intent(in) :: t(n),c(n),x(m)
real(RKIND), intent(out) :: y(m)
real(RKIND), intent(inout) :: wrk(n)
! ..local scalars..
integer :: i,j,k1,k2,k3,l,l1,l2,nk1,nk2,kk
real(RKIND) :: ak,arg,fac,tb,te
! ..local arrays ..
real(RKIND) :: h(MAX_ORDER+1)
logical :: nonflat
! before starting computations a data check is made. if the input data
! are invalid control is immediately repassed to the calling program.
ier = FITPACK_INPUT_ERROR
if (nu<0 .or. nu>k) return
if (m<1) return
kk = k-nu
ier = FITPACK_OK
! fetch tb and te, the boundaries of the approximation interval.
k1 = k+1
k3 = k1+1
nk1 = n-k1
tb = t(k1)
te = t(nk1+1)
! the derivative of order nu of a spline of degree k is a spline of
! degree k-nu,the b-spline coefficients wrk(i) of which can be found
! using the recurrence scheme of de boor.
l = 1
wrk(1:nk1) = c(1:nk1)
if (nu/=0) then
nk2 = nk1
de_boor: do j=1,nu
ak = k1-j
nk2 = nk2-1
l1 = l
do i=1,nk2
l1 = l1+1
l2 = l1+k1-j
fac = t(l2)-t(l1)
if (fac>zero) wrk(i) = ak*(wrk(i+1)-wrk(i))/fac
end do
l = l+1
end do de_boor
endif
nonflat = nu==0 .or. k/=nu
l = k1
l1 = l+1
k2 = k1-nu
j = 1
! main loop for the different points.
user_points: do i=1,m
! fetch a new x-value arg.
arg = x(i)
! check if arg is in the support
if (arg < tb .or. arg > te) then
select case (e)
case (OUTSIDE_EXTRAPOLATE)
! continue like any other point
case (OUTSIDE_ZERO)
y(i) = zero
cycle user_points
case (OUTSIDE_NOT_ALLOWED)
ier = FITPACK_INSUFFICIENT_STORAGE
return
end select
endif
! search for knot interval t(l) <= arg < t(l+1)
do while (.not.(arg>=t(l) .or. l1==k3))
l1 = l
l = l-1
j = j-1
end do
! ++
do while (.not.(arg<t(l1) .or. l==nk1))
l = l1
l1 = l+1
j = j+1
end do
if (nonflat) then
! evaluate the non-zero b-splines of degree k-nu at arg.
h = fpbspl(t,n,kk,arg,l)
! find the value of the derivative at x=arg.
y(i) = dot_product(h(1:k2),wrk(l-k:l-nu))
else
! if nu=k the derivative is a piecewise constant function
y(i) = wrk(j)
endif
end do user_points
return
end subroutine splder